Initialization

Parameters

Cp<-1005 #specific heat capacity
R=287.058 # specific gas constant for dry air (J/(kg·K))
lambda=4
wd<-"C:/Users/marie/OneDrive/Documenten/uhi_max/"

svf<-stack(paste0(wd,"inst/Grids_veg_svf/svf_utrecht_1m.grd"))
fveg<-stack(paste0(wd,"/inst/Grids_veg_svf/greenness_utrecht_smooth_500m.grd"))
fveg<-projectRaster(fveg,crs=crs(svf))

IUtrecht23<-data.frame(lon=52.079,lat=5.139)
coordinates(IUtrecht23)<- ~lat+lon
crs(IUtrecht23)<-crs("+init=epsg:4326")

mapview(fveg) + mapview(IUtrecht23)
## Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ; 
## the supplied Raster* has 1348088 
##  ... decreasing Raster* resolution to 5e+05 pixels
##  to view full resolution set 'maxpixels =  1348088 '
IUtrecht23_fveg<-raster::extract(fveg,IUtrecht23)
## Warning in .local(x, y, ...): Transforming SpatialPoints to the CRS of the
## Raster
IUtrecht23_svf<-raster::extract(svf,IUtrecht23)
## Warning in .local(x, y, ...): Transforming SpatialPoints to the CRS of the
## Raster

From 10min data to Daily UHImax parameters

From the autmaomatic weather stations (AWS) in the rural area the following data is required:

The 10 minute data from global radiation, wind speed, temperature and pressure is used to calculate the UHImax. The hourly relative humidity, precipitation and mean wind are used to filter out synoptic situations with frontal system and fog. Under these different conditions there is no clear relation between the AWS and city temperature.

cabauw_10min<-fread(paste0(wd,"inst/Cabauw_meteo/Cabauw_T_S_W_2016_201809.csv"))
cabauw_10min$IT_DATETIME<-as.POSIXct(cabauw_10min$IT_DATETIME,format="%Y%m%d_%H%M00_000000")

Cabauw.S<-subset(cabauw_10min[which(cabauw_10min$DS_CODE=="348_S_a"),],
                 select = c("IT_DATETIME","DS_CODE","TOS.Q_GLOB_10"))
Cabauw.T<-subset(cabauw_10min[which(cabauw_10min$DS_CODE=="348_T_a"),],
                 select = c("IT_DATETIME","DS_CODE","TOT.T_DRYB_10"))
Cabauw.W<-subset(cabauw_10min[which(cabauw_10min$DS_CODE=="348_W_a"),],
                 select = c("IT_DATETIME","DS_CODE","TOW.FF_SENSOR_10"))

S<-calc_S(Cabauw.S$IT_DATETIME,Cabauw.S$TOS.Q_GLOB_10)
DTR<-calc_DTR(Cabauw.T$IT_DATETIME,Cabauw.T$TOT.T_DRYB_10)
U<-calc_U(Cabauw.W$IT_DATETIME,Cabauw.W$TOW.FF_SENSOR_10)

Meteo_params<-data.frame(DTR,"S"=S$S[1:length(DTR$start)],"U"=U$W)
Meteo_params$start<-as.POSIXct(Meteo_params$start)
Meteo_params$stop<-as.POSIXct(Meteo_params$stop)
saveRDS(Meteo_params,paste0(wd,"inst/Rdata/Meteo_params_cabauw.rds"))

Data preprocessing

Within the city of Utrecht we have carefully selected three stations from the Wunderground network:

These stations have been measuring the city temperature for the past couple of years. The data can be downloaded using download_time_seq, note that the temperature is in Farenheid and not degrees Celcius. The data has been filtered and interpolated to 10minute timestamps. Unrealistic high and low temperatures are set to NA. In case more than 8 identical measurements in a row occur the value is also set to NA.

cabauw_parms<-readRDS(paste0(wd,"inst/Rdata/Meteo_params_cabauw.rds"))

cabauw_P<-fread(paste0(wd,"inst/Cabauw_meteo/Cabauw_P_10min.csv"))
cabauw_P$IT_DATETIME<-as.POSIXct(cabauw_P$IT_DATETIME,format="%Y%m%d_%H%M00_000000")

cabauw_T<-fread(paste0(wd,"inst/Cabauw_meteo/Cabauw_T_10min.csv"))
cabauw_T$IT_DATETIME<-as.POSIXct(cabauw_T$IT_DATETIME,format="%Y%m%d_%H%M00_000000")
cabauw_PT<-merge(cabauw_P,cabauw_T,by="IT_DATETIME")

air_density<-calc_U(cabauw_PT$IT_DATETIME,(cabauw_PT$TOA.P_NAP_MSL_10*100)/(R*(cabauw_PT$TOT.T_DRYB_10+273.15)))
names(air_density)<-c("start","stop","rho")
cabauw_parms<-merge(cabauw_parms,air_density,by=c("start","stop"))
cabauw_parms$S_new<-cabauw_parms$S/(cabauw_parms$rho*Cp)
cabauw_parms$meteo<-(((cabauw_parms$S_new)*cabauw_parms$DTR^3)/cabauw_parms$U)^(1/lambda)

# svf_grn<-readRDS(paste0(wd,"inst/Rdata/wunderground_svf_grn.rds"))
# df.svf_grn<-data.frame(svf_grn)
# df.svf_grn<-df.svf_grn[df.svf_grn$Station.ID %in% c("IUTRECHT196","IUTRECHT376","IUTRECHT299"),]

Calculate UHI max parameters

Setting all parameters to required to calculate the UHI with the formula of Theeuwes (2017).

# for(i in 1:length(df.svf_grn$Station.ID)){
# STN<-df.svf_grn$Station.ID[i]
STN<-"IUTRECHT23"
stn<-readRDS(paste0(wd,"inst/Wunderground/Filtered/IUTRECHT23_filtered.rds"))
wur_cabauw<-merge(x=stn$interpolated_time,y=cabauw_T,by.x="new_time",by.y="IT_DATETIME")
wur_cabauw<-wur_cabauw[complete.cases(wur_cabauw),]
UHI<-calc_UHImax(time = wur_cabauw$new_time,
                        Tcity = wur_cabauw$T_int,
                        Tref = wur_cabauw$TOT.T_DRYB_10)
UHI$Tcity<-calc_U(time=wur_cabauw$new_time,wur_cabauw$T_int)$W
UHI$Tref<-calc_U(time=wur_cabauw$new_time,wur_cabauw$TOT.T_DRYB_10)$W
# UHIcalc_stn<-cbind(cabauw_parms,UHIcalc)
UHI<-merge(UHI,cabauw_parms,by=c("start","stop"))
UHI$svf<-as.numeric(IUtrecht23_svf)
UHI$fveg<-as.numeric(IUtrecht23_fveg)
UHI$Cp<-Cp
UHI$stn<-STN


# p<-ggplot(UHI,aes(UHImeasured,meteo))+geom_point()+geom_abline()
# ggsave(p,filename=paste0(wd,"inst/UHImax/fig/",
#                 STN,".png"))
write.table(UHI,paste0(wd,"inst/UHImax/",
                       STN,"_UHIparams.txt"),
            row.names = FALSE,
            col.names = TRUE,
            sep=",")
# }

Subset based on weather conditions

####################Select control days
meteo.df<-fread(paste0(wd,"inst/Cabauw_meteo/rain_wind_rh_hour.csv"))
meteo.df$IT_DATETIME<-as.POSIXct(meteo.df$IT_DATETIME,format="%Y%m%d_%H%M00_000000")
time<-meteo.df$IT_DATETIME
rh<-meteo.df$BGH.U
rain<-meteo.df$BGH.Q_RH
wind<-meteo.df$BGH.FH

days_subset<-uhi_sub(time=time,wind=wind,rain=rain,rh=rh)
######################
uhimax_files<-list.files(paste0(wd,"inst/UHImax/"),
                         full.names = TRUE,pattern=".txt")
uhimax_files<-lapply(uhimax_files,fread)
uhimax_Utrecht<-do.call("rbind",uhimax_files)

uhimax_Utrecht<-merge(days_subset,uhimax_Utrecht,by=c("start","stop"))
# uhimax_Utrecht<-uhimax_Utrecht[which(uhimax_Utrecht$Rain==TRUE),]
# uhimax_Utrecht<-uhimax_Utrecht[which(uhimax_Utrecht$Wind==TRUE),]
# uhimax_Utrecht<-uhimax_Utrecht[which(uhimax_Utrecht$rh==TRUE),]
uhimax_Utrecht<-uhimax_Utrecht[which(uhimax_Utrecht$Select==TRUE),]
uhimax_Utrecht<-uhimax_Utrecht[which(uhimax_Utrecht$Tref>17),]
uhimax_Utrecht<-uhimax_Utrecht[which(uhimax_Utrecht$stn %in% c("IUTRECHT23")),]
# uhimax_Utrecht<-uhimax_Utrecht[which(uhimax_Utrecht$stn %in% c("IUTRECHT23","IUTRECHT196","IUTRECHT376","IUTRECHT299")),]

Correlations with the three meteo params

uhi_u<-ggplot(uhimax_Utrecht,aes(U,UHImeasured,colour=DTR)) +geom_point() + scale_color_gradientn(colours = terrain.colors(20))
uhi_DTR<-ggplot(uhimax_Utrecht,aes(DTR,UHImeasured,colour=S_new)) + geom_point()+scale_color_gradientn(colours = heat.colors(20))
uhi_S<-ggplot(uhimax_Utrecht,aes(S_new,UHImeasured,colour=U)) +geom_point() + geom_point()+scale_color_gradientn(colours = topo.colors(20))
ggsave(uhi_u,filename = paste0(wd,"inst/UHImax/fig/uhi_u.png"))
## Saving 7 x 5 in image
ggsave(uhi_DTR,filename = paste0(wd,"inst/UHImax/fig/uhi_DTR.png"))
## Saving 7 x 5 in image
ggsave(uhi_S,filename = paste0(wd,"inst/UHImax/fig/uhi_S.png"))
## Saving 7 x 5 in image
uhi_u

uhi_DTR

uhi_S

ggplot(uhimax_Utrecht,aes(UHImeasured,(2-svf-fveg)*meteo,colour=factor(stn))) +geom_point() +geom_abline() +xlim(0,10)+ylim(0,10)

caret::R2((2-uhimax_Utrecht$svf-uhimax_Utrecht$fveg)*uhimax_Utrecht$meteo,uhimax_Utrecht$UHImeasured)
## [1] 0.2123518
caret::RMSE((2-uhimax_Utrecht$svf-uhimax_Utrecht$fveg)*uhimax_Utrecht$meteo,uhimax_Utrecht$UHImeasured)
## [1] 0.9565242

Using the formula of Theeuwes for Utrecht

The formula of Theeuwes(2017) was derived for:

Outside of this range the formula has not been tested.

svf.r<-resample(svf,fveg,method="bilinear")
values(svf.r)[values(svf.r)<0.2 | values(svf.r)>0.95] = NA
values(fveg)[values(fveg)<0 | values(fveg)>0.4] = NA

city_center<-extent(c(5.1,5.13,52.08,52.098))
city_center<-raster(city_center)
values(city_center)<-1
crs(city_center)<-CRS("+proj=longlat +ellps=WGS84 +datum=WGS84 +no_defs")
city_center<-projectRaster(city_center,crs=crs(fveg))

st<-(2-svf.r-fveg)*mean(uhimax_Utrecht$meteo)
st<-crop(st,city_center)
mapview(st)